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Abstract 

A real-space representation of the current response of many-electron systems with possible ap- 
plications to x-ray nonlinear spectroscopy and magnetic susceptibilities is developed. Closed ex- 
pressions for the linear, quadratic and third-order response functions are derived by solving the 
adiabatic Time Dependent Current Density Functional (TDCDFT) equations for the single-electron 
density matrix in Liouville space. 
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I. INTRODUCTION 



Time Dependent Current Density Functional Theory (TDCDFT) offers a computation- 
ally tractable framework for computing currents and response functions of many-electron 
systems in response to external electric and magnetic perturbations- The time-dependent 
linear paramagnetic susceptibilities are calculated as the linear response of a noninteracting 
Kohn-Sham reference system to an effective vector-potential, which consists of the external 
field, together with the Hartree and the exchange-correlation contributions 1 . For the sake 
of computing current-related properties, it is natural to consider the effective potential to 
be a functional of current and charge density rather than charge density alone, as in stan- 
dard Time Dependent Density Functional Theory (TDDFT) 2 -. Another reason for applying 
TDCDFT to crystals is connected with the recent argument that there is a one-to-one cor- 
respondence between time-varying periodic potentials and the current density but not with 
the charge density 3 . TDCDFT was successfully used for calculating the polarizability of 
conjugated polymers^. Current density functional theory (CDFT) yields exact response 
functions^ to static external potentials and TDCDFT is thus expected to provide reasonable 
approximations for time-dependent current properties. 

The linear magnetic susceptibility (the response of the current to an external vector- 
potential) of the Kohn-Sham non-interacting system has been calculated using the local 
current density exchange correlation kernel for the electron gas 1 -. However, computing the 
response functions of the interacting system (where the quasiparticle energies cannot be 
expressed as differences between Kohn-Sham energy levels), requires the solution of a chain 
of integral equations 1 , whose computational cost rapidly increases with the nonlinear order 
of the response. In this letter we compute current response functions by recasting the TD- 
CDDT equations in terms of the reduced single electron density matrix for an N electron 
system p(r, r 1; t) = J2n=i VVi( r j t)tp^(r 1 , t) 7 , where ip n ( r ) are the Kohn-Sham orbitals. Closed 
expressions are derived for the linear, quadratic and third-order response function o 2 ^ 1 ? (the 
response of the total polarization current to the external electric field) by solving an eigen- 
value equation in Liouville space. The quasiparticle frequencies are not given simply as 
differences of Kohn-Sham orbital energies^. Current response functions should also be par- 
ticularly suitable for computing the resonant nonlinear response to x-ray fields 10 . Due to 
the short wavelength, the dipole approximation does not generally hold, and x-ray suscepti- 
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bilities may be expressed in terms of multitime correlation functions of currents and charge 
densitie a 10 ! 11 . TDCDFT thus provides a natural direct computational approach for nonlinear 
x-ray spectroscopy. 



II. TIME DEPENDENT KOHN-SHAM CURRENT-DENSITY FUNCTIONAL 
EQUATIONS FOR THE SINGLE ELECTRON DENSITY MATRIX 

TDCDFT maps the original system of interacting electrons onto an effective system of 
non-interacting electrons subjected to an exchange-correlation scalar and vector potentials^, 
constructed to yield the same charge density and current profiles as the interacting system. 
The Kohn-Sham TDCDFT equations of motion for the time-dependent density matrix are 

■ dp(r,ri,t) 

1 dt 



H° KS (r, t) - H% s (n, f) p(r, r 1; t) - (j(r, t) A(r, t) - j*(r 1; t)A*(r 1: *)) 



+ (^ M)nM) _^ (ri)tK(ri)t) ) ; (1) 

where e(m) is the electron charge (mass), and we set h — 1. The time- dependent charge 
density and paramagnetic electronic current are given by 

n(r, t) = p(r, r, t); j(r, t) = -— [(V ri - V r )5p(r, n, t)] r=Ii . (2) 

The observed (physical) current, which enters the continuity equation, is given by 

e 2 

J(r, t) = j(r, t) - 2^A(r, t)n(r, t); (3) 

Hft S is the Kohn-Sham Hamiltonian H^ s (r,t) = H^ s [n{r, t), j(r, t)]. The remaining terms 
in Eq. (|1J) represent the coupling to an external vector-potential. 

H° KS [n(r,t),j(r,t)} = -^V 2 + U KS (r, t), (4) 

where the potential Uks(. t , t) = UKs[ n ( r , t), j(r, t)} is 

U KS (r, t) = -^V r A xc (r, t) - ^A xc (r, t)V T + U° KS (r, t) + U (r). (5) 

The exchange-correlation vector potential A xc [n(r,t), j(r,t)] and the Kohn-Sham scalar 
external potential U^ s [n(r,t), j(r,t)] are functionals of both the charge density and the 
current 1 . The scalar potential is given by: 

U° KS (r,t) = f dn " (ri>t) f + U xc (r,t); (6) 



Uq(t) is the field created by nuclei and U xc [n(r,t),j(r,i)](r) is the exchange-correlation po- 
tential in the adiabatic approximation. The time-dependent external potential is U ex t(r, t) = 
U (r) at time t < t and U ext (r,t) = U (r) + Ui(r, t) for t > t . A xc adds a magnetic field 
induced by the exchange-correlation interaction between electrons. Note that unlike the 
paramagnetic canonical current, J is gauge invariant 1 . 

The stationary solution of Eq. (JI} gives the ground state single electron density matrix 
p(r, ri) which carries no current. We then set p(r, rx,t) = p(r, ri) + Sp(r, rx, t) where 
5p represents the changes induced by Ui(r,t). Its diagonal elements 8n(r) = 5p(r,r,t) 
give the changes in charge density, whereas the off-diagonal elements represent the changes 
in electronic coherences between two points. The physical current may be obtained by 
expanding Sp in powers of the external vector potential A(r, t): 5p = 8p\ + 5p 2 + ... and 
solving Eqs. (JIJ and (J2|) self-consistently for 5p order by order. To that end we first recast 
Ujcs[ n ( r :t),j(r,t)] as a functional of the paramagnetic current j(r, t) alone. This is done 
by substituting the total current J(r, t) (Eq. (jHJ)) into the continuity relation between the 
charge density and the total current 

8n{r, t) = - X - J* V r J(r, r)dr. (7) 

Solving Eq. Q self-consistently for the charge density n(r, t) in terms of the paramagnetic 
current j(r, t), and substituting it in Uks[ii(t, t),j(r, t)] eliminates the explicit dependence 
on the charge density. Expanding Uks around p to second order in j, we obtain 

U KS [i(r,t)](r) = U° KS [p](r) + J dv x {j^^ + flM*, r i)J J( r i> *) 
+ / dr 1 J dr 2 g~ xc [p\(r,v 1 ,r 2 )j(r 1 ,t)j(v 2 ,t); (8) 

where f xc \p)(r, r x ) and g~ xc [p](r, r 1; r 2 ) are the first and the second order adiabatic exchange- 
correlation kernels. We have made the commonly used adiabatic approximation where we 
assume that the kernels are time-independent. 

III. QUASIPARTICLE REPRESENTATION OF X-RAY NONLINEAR RE- 
SPONSE FUNCTIONS 

We next separate 8p into an electron-hole, interband, (£) and intraband (T) components 
5p(r, Tx,t) = £(r, ri,£) + T(£(r, r 1; t)). It follows from the idempotent property of p, that 
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T is uniquely determined by The matrix elements of £, unlike those of 5p, constitute 
independent variables, that can be used to construct a quasiparticle representation. 

The quasiparticle spectrum is obtained by solving the linearized Kohn-Sham eigenvalue 
equation 

L£ a (r, r x ) = Q a £ a (r, r ± ), (9) 

where 

L£ a (r,ri) = -— — + P( r > r i) (f dr i J dr 2 ^Jp](r, ri, r 2 )£ a (r 1; r 2 )) 

- p*(r u r) (j dr J dr a /£|fl(r Xj r, r 2 )£ a (r 2 , r)) , (10) 

fL[p}( r , r i ; r 2) = /zc [p](r,n)5(ri - r 2 )w(ri,r 2 ), (11) 
and -0 is the velocity operator 

«(r,n) = ~ (Vn-V r ). (12) 

The quasiparticle modes £ Q come in pairs a = ±1,±2, ... with = Q a . Their or- 
thonormality and algebraic properties make it possible to expand an arbitrary interband 
matrix^ in the form 

£(r,ri,t) = X;£ Q (r,ri)z a (t), (13) 

a 

where z a (t) = where the scalar product of any two interband matrices £ and r\ is 

defined by^ 

<ei»7> =J dr J dr'p[^,v](r,r')5{r-r'). (14) 

The bra (ket) notation underscores the similarity with Dirac's Hilbert space notation. 
z -a(t) = z a{t) constitute complex oscillator amplitudes. We shall denote their perturbative 
expansion in the external vector-potential A(r, t) by z^\ 

The paramagnetic current is calculated by taking the expectation value of the velocity 
v(r, ri) with respect to the time-dependent density matrix 

a af3 

+ \ ^}c t ,M( r ^') z ^{ t )zp{t)z 1 {t), (15) 

6 a/37 



where (a, (3, 7 = ±1, ±2, . . .) and we only retained terms that contribute to the second order 
response. Here 



ic f 

j«(r) = J dr x [5(t - n)(V ri - V r )e a (r, n)] 



(16) 



j«/?(r) = "^/ rfr i 5(r-n)(V ri - V r ) JdT 2 {{S{T 2 -T) 

- 2n(r,r 2 )) / dr 3 (£ a (r 2 , r 3 )^(r 3 , r x ) + ^(r 2 , r 3 )£ Q (r 3 , n)) 

ja/3 7 ( r ) = ^J dri ^( r - r i)(V ri - V r ) y dr 2 (£ Q (r,r 2 ) 
/ dr 3 (^ a (r2,r 3 )^(r3,ri) + ^(r 2 , r 3 )£ Q (r 3 , n)) 



(17) 



(18) 



The collective electronic oscillator (CEO) expansion for the charge density n(r,t) is given 
by 9 

Sn(r,t) = $^n a (r)z a (t) + )-^n aP (r)z a (t)zp(t) 

a 1 a/3 

1 



d a/3 7 

a,/3, 7 = ±l,±2,... 



(19) 



The coefficients of this expansion are given by Eqs. (29)- (31) in RefA 

The equations of motion for £ can be obtained from Eq. by expressing the density 
matrix in terms of £ and T(£) 9 . Equations of motions for z^'(t) are derived in terms of 
A(r,i) and j(r) by substituting the mode expansion of £ into these equations. Substituting 
the solution z^\t) into Eq.([T5)). gives 

ft 



t)= I drfdrt J2 x£i.(*.'",r,ri).4 Al (r 1 ,T), 



(20) 

fi=x,y,z 

where x is the linear paramagnetic susceptibility, and A s , Ai are Cartesian tensor compo- 
nents. 

We further introduce the observed susceptibility x \1 (t > r > r i r i) defined by replacing the 
paramagnetic current j(r, t) with the physical current J(r, t) in Eq. (J2DJ)- Substituting 
j (1) (r,i) from Eq. © into Eq. JUJ), we obtain 



X^^.r.ri) = Xv(^,^r,ri) - ^-7i(ri)<J(ri - r)5(r - t)5 x ». 
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(21) 



We next derive equations of motions for z^{t) in terms of A(r, t) and j a (r), by substituting 
the CEO expansion of £ from Eq. (fT3*j) into the equation of motion for £, which can be 
obtained from Eq. ((TJ by expansion of the density matrix on £ and T(£) 9 . Solving these 
equations we substitute z^(t) into Eq.(fT5]). 

«2>(<*,r,r0 = £ 2 "^ (r) ^ (r,) . (22) 

o=l,2,... 

Optical with x-ray signals are most conveniently expressed using response functions which 
connect the polarization with the electric field. For example, the linear response function 
a \iX B (A r ' r ' r i) ^° fi rs ^ or der in the external field E(r, t) is defined as: 

Pg ) (r,t)=/* drfdr, £ a^(t, r, r, r), (23) 

where P(r, i) is the total polarization^. a^\t, r, r, ri) can be obtained from Eq. (J20|) by 
noting that j is connected to J through Eq. (|3*|): J is connected to P through^ P(r, t) = 
f^oo drJ(r, t) and A(r, u;) = — zcE(r, uo)/uo. Using these relations we obtain^ 

•l-n 

a {n) {uj, r,r n , . . . , n, w n , . . . , u x ) = x (n) (^,r,r n , . . . ,Ti,u n , . . (24) 

cji^ . . . k; n u; 

this gives for the linear response 

(!) 1 \ 1 (!) / \ 

r ' ri ) = 7aX\i\,( w = w i' r ' ri ) 

to 1 

1, / r 2ti^-(r)j^( ri ) e 2 \ 

= t Ji,... ai-^ 2^' !(ri)<5(ri " r)<5 "' j ' (25) 

where the first term in the bracket is Xx Xi( u i r > r i); an< ^ ^( r ) * s the g roun d state charge 
density. Eq. (J2*5j) provides a microscopic algorithm for computing the Kubo formula^; all 
quantities are obtained from the quasiparticle modes. 
To calculate the second order response function 

PxJ(r,t) = \ I dn j dr 2 I drt f dr 2 Ex l (ri,T 1 )E\ 2 (r 2 ,T 2 ) 

Z J— 00 J— 00 J J 

^A?A a A.(*» 7 i J r 2,r,r a ,r 2 ), (26) 

we introduce the second-order exchange-correlation kernel g xc , obtained by expanding the 
exchange correlation potential UKg\j(r,t)](r) to the second order by <5j (Eq. (jHJ)) 

g xc [n]{r, r 1; r 2 , r' 2 , r 4 ) = g~ xc [n](r, r a , r 3 )5(ri - r 2 )5(r' 2 - r 4 ){)(r 1 , r 2 )v(r / 2 , r 4 ). (27) 



Repeating the procedure used for to the next order we obtain the second order param- 
agnetic susceptibility 

~(2) , X o V- ^(-«/37)( r ' r l' r 2)ia S ( r )j-M ri )A( r 2) S a S /3 



«A s /37 



(Q a — Ui — UJ 2 )(^i3 — wi)(0 7 — U0 2 ) 



+ y 3t' a p(T)jt(Tl)j%(T2)s a S | (r) J^ 1 (tQ (r 2 ) S a Sp 

+ Z^T^i vFS V' a 'A7 = ±1,±2,..., 28 

where s a = sign(a). V^(_ aj g 7 ) is obtained by substituting the exchange-correlation kernels 
f xc from Eq. (jTTjl and g xc into the expression for V g ^_ af3 ^ given in Ref.— . 

Similar to the linear response function, the second-order response function is finally ob- 
tained by expanding of the charge density n(r, t) in the modes (Eq. ((IS)) ): 

e 2 2^(r)n a (n) 



2mc \ ^ Vt 2 a - (uji + uj 2 ) 2 

(29) 



a fi a - (^1 + 



Higher response functions can be computed similarly^. The third-order response function is 
given in Appendix [X] 



IV. DISCUSSION 



To get the high-order paramagnetic susceptibilities in the standard Hilbert space TD- 
DCDT approach one needs to solve self-consistently a chain of integral equations for each 
order 1 . The linear paramagnetic susceptibility in the standard Hilbert space TDDCDT ap- 
proach is given by Eqs. (8)- (9) in Ref.—. In contrast, the closed expressions for the linear 
(Eq. (J23)), second-order (Eq. (J2EJ)) and the third-order (Eq. (jA12)l ) susceptibilities derived 
in this paper use the CEO representation in Liouville space. 

Correlation-function expressions for the linear, second-order and third-order x-ray re- 
sponse functions were derived in Eqs.(Bl),(B2),(B3a)- (B3d) in RefA Eqs. (J23J) 
and (jAl|) express these TDCDFT response functions in the CEO representation, and pro- 
vide a computational scheme for nonlinear x-ray response functions. 
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TDDFT exchange-correlation functional are better developed and more widely used 
than their TDCDFT counterparts. TDDFT currents can be obtained by simply modifying 
Eqs. ff llfl by setting v = 1, and using the TDDFT exchange-correlation kernels^ where the 
scalar exchange-correlation potential depends only on charge density. 

Finally we note that this work can be extended to include non-adiabatic exchange- 
correlation potentials, as outlined recently for the linear response^. In general, the exchange- 
correlation potential and exchange-correlation kernels are time-dependent 14 . This time de- 
pendence has been neglected within the adiabatic approximation used here. If we relax this 
approximation, the eigenvalue equation for the Liouville superoperator L, Eq. (jHJ), should 
be replaced by^ 

L(Q a )Ur,r') = Q a Ur,r'). (30) 

Methods for solving Eq.(j20J) using the frequency- dependent functional of Gross and Kohni^ 
were described in^- 



Acknowledgements 

This article is based upon work supported by the National Science Foundation under 
grant number (CHE-0132571). This support is gratefully acknowledged. 



APPENDIX A: THE THIRD-ORDER RESPONSE 



For the third-order response function we obtain 



X ( 2x 2 \ 3 \ s w 2 , u 3 , r, n, r 2 , r 3 ) 



1 



2tt 



: — (F[j^(r)n(r 3 )j Al (ri)]5(r 2 - r 3 )5 A2 A 3 + F[j^(r)n(r 3 )j Al (ri)]5(r 2 - r,)^ 



2mc 

+ F[n(r)j A2 (r 3 )j Al (r, )}Mv - r,,Rw V b 

/ .2 \ 2 



+ 



2mc 



i? (1) (wi + u 2 + uj 3 , r, r 2 )5(r - r 3 )5(r 2 - r 1 )5x 3 \J\ 2 \ 1 



(Al) 



where R^>(u>i + cj 2 + to 3 ,r, r 2 ) is the linear density-density response given by Eq. (46) in 
Ref. 9 ; xi 3 1 ) A 2 A 3 A s (o; 1 ,a; 2 ,a; 3) r,r 1) r 2) r 3 ) is the third-order paramagnetic susceptibility given 



by the r.h.s. of Eqs. (C23)- (C31) in Ref.— by replacing p by j: 

perm 

X ( xLx 3 xM^2,^,r,r'yy')= £ {x? ) +X?l+X?l I + ---xPin), (A2) 



where 



-(3) = y J^(r)j% 7 (rQ^(^)j^(r w )s ag/3g7 

-(3) = y j^(r)y g( _^ ) j^(r / )j^(r w )j^(r w )s a s /3 s 7 s, 

^ = E tf; — ; TTTik ; TXTci rr> ( A5 ) 



^ (fi Q -ui-ui — u; 3 )(fis - a> 2 - ^ 3 )(^ 7 - w 3 ) 
Q/ 3 7 5 (^a - ^1 - ^2 - ^ 3 ) - ^i) (fi 7 - uj 2 - uj 3 ) (Q 5 - V3) ' 



-(3) \- 2^ ( _ 0/37 )j^(r)j^(r / )j^(^)j^(r w ) gaS/3g7S5 

= E 775 ; ; — —\7T, Tvf; ; txtt, rr> ( A6 ) 



~(3) 

Xv = 

y 2^ ( _^ 7 )^(_ 7 ^)j^(r)j^(rQj^Kr // )j^( rW ) g ^^7^^ (A?) 

^ c^5 (^a ~~ ^1 ~~ ^2 - ^3)(^/3 - ^l)(^ 7 - ^2) (^5 - ^3)' 1 

-(3) _ ^ ^(r)j^(r , )j^(r ,/ )A(r //, )aaa/ ? a 7 / » Q \ 

Xyn — 7r5 TTT?) i , w n Ti l Ay J 



a/37 



(fi a - Wi)(tyg - W 2 - ^3) (^7 - ^3) 



-(3) v ^(^^(-^jJ-UrQj^lr^J^K^)^^^^ , Aim 

XvHI ^ 5 (n a - u,i)(fy - u 2 - u 3 )(n 7 - u 2 )(n s - u 3 y 1 J 

Here v = a, (3, 7, 5, 77 = ±1, ±2, . . . and fi^ is positive (negative) for all ^ > (z/ < 0) 
according to the convention Q^ u = —Q u . 

F[j As ( r ) n ( r 3)j Al ( r i)] i s determined by the r.h.s. of Eq. ([29)1 by replacing j Al (i"i) by n(r 3 ): 

F[j A Xr)n(r 3 )j Al (ri)] = - ^'^ tt [^J^.^.r,^^) 

' ' (AH) 



^ 2fi Q n a (r)n(r 3 ) 



where 

(2) , , ^(-a/3 7 )( r ; n. r3)j As (r)^-/3(r3)j-7(ri)g a 5/3 
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+ l^Tn vo V 7 = ±1, ±2, ■ • • , (A12) 

where s a = sign(a). 
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